{
"cells": [
{
"cell_type": "markdown",
"id": "3d4fc10a",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"# Error-tolerant BS-based circuit"
]
},
{
"cell_type": "markdown",
"id": "8cdc5601",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"In this notebook, we aim at presenting an alternative form of \"generic\" interferometers using a BS-based structure as a building block rather than a MZI (Mach-Zehnder interferometer) as a building block. BS-based circuits seem to be more tolerant to manufacturing errors and losses according to Fldzhyan et al. *Optimal design of error-tolerant reprogrammable multiport interferometers*, [Optics Letters](https://doi.org/10.1364/OL.385433), 45(9):2632–2635 (2020)"
]
},
{
"cell_type": "markdown",
"id": "561b1415",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"## Introduction "
]
},
{
"cell_type": "markdown",
"id": "94227072",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"### Parameters"
]
},
{
"cell_type": "markdown",
"id": "897a99d7",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"The goal is to successfully implement a random unitary $U_{target}$ by varying parameters of a fix \"generic\" interferometer $U_{interf}$. For instance, the Reck decomposition [2] has MZI as building block and by varying the angle of variable phase shifts, any unitary can be implemented. However, as we will see, this is only true when the MZI building block have perfect 50:50 beam-splitters which is not realistic in practice due to manufacturing errors. To see this, we will introduce a parameter $\\alpha$ which characterise the error caused by the imbalance of the static BS due to imperfect realisation. And we will investigate how this impacts the correct implementation of any unitary (drawn from the Haar measure) by minimising the infedility (or maximising the fidelity) between the target unitary $U_{target}$ and its implementation $U_{interf}$.\n",
"We will compare this with the BS-based approach [1]."
]
},
{
"cell_type": "markdown",
"id": "5f57bc2b",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"#### Initialisation"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "190e653d",
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"import datetime\n",
"import time\n",
"import numpy as np\n",
"from scipy.optimize import basinhopping\n",
"import random\n",
"import perceval as pcvl\n",
"import perceval.lib.phys as phys\n",
"import os\n",
"import math\n",
"from tqdm.notebook import tqdm_notebook"
]
},
{
"cell_type": "markdown",
"id": "bcd40037",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"We start by defining all parameters used in the program later on:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "09dfc317",
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"N=5\n",
"n_try=10\n",
"n_iter=3\n",
"n_process=8\n",
"angle_min=-15\n",
"angle_max=30\n",
"angle_step=2\n",
"n_unitary=300\n",
"logfilebs='bsbasednotebook-opt'\n",
"logfilemzi='mzibasednotebook-opt'"
]
},
{
"cell_type": "markdown",
"id": "98d2cc35",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"## 1. Perceval implementation of the BS-based interferometer"
]
},
{
"cell_type": "markdown",
"id": "e6c9c678",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"We start by describing the circuit as defined in [1] using `pcvl.Circuit.generic_interferometer()`. The circuit is built with an arrangement of single static beam-splitters with a variable phase shifter on the same leg. It also starts with phase-shifters on each mode at the beginning."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "0081f2c8",
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"text/plain": "<IPython.core.display.HTML object>",
"text/html": "<svg width='1412.000000' height='300.000000' viewBox='-6.000000 0.000000 1406.000000 300.000000'><path d=\"M 10 25 L 25 25\" fill=\"none\" stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\"/>\n<text x=\"0.000000\" y=\"25.000000\" font-size=\"12.000000\" text-anchor=\"start\">0</text>\n<path d=\"M 10 75 L 25 75\" fill=\"none\" stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\"/>\n<text x=\"0.000000\" y=\"75.000000\" font-size=\"12.000000\" text-anchor=\"start\">1</text>\n<path d=\"M 10 125 L 25 125\" fill=\"none\" stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\"/>\n<text x=\"0.000000\" y=\"125.000000\" font-size=\"12.000000\" text-anchor=\"start\">2</text>\n<path d=\"M 10 175 L 25 175\" fill=\"none\" stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\"/>\n<text x=\"0.000000\" y=\"175.000000\" font-size=\"12.000000\" text-anchor=\"start\">3</text>\n<path d=\"M 10 225 L 25 225\" fill=\"none\" stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\"/>\n<text x=\"0.000000\" y=\"225.000000\" font-size=\"12.000000\" text-anchor=\"start\">4</text>\n<polyline points=\"25 25 75 25\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"30 40 39 40 53 10 44 10 30 40 39 40 30 40\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"47.000000\" y=\"38.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=3.223159</text>\n<polyline points=\"25 75 75 75\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"30 90 39 90 53 60 44 60 30 90 39 90 30 90\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"47.000000\" y=\"88.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=4.428628</text>\n<polyline points=\"25 125 75 125\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"30 140 39 140 53 110 44 110 30 140 39 140 30 140\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"47.000000\" y=\"138.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=5.440784</text>\n<polyline points=\"25 175 75 175\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"30 190 39 190 53 160 44 160 30 190 39 190 30 190\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"47.000000\" y=\"188.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=2.466231</text>\n<polyline points=\"25 225 75 225\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"30 240 39 240 53 210 44 210 30 240 39 240 30 240\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"47.000000\" y=\"238.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=3.407489</text>\n<polyline points=\"75 25 103 25 122 44\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"128 44 147 25 175 25\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"75 75 103 75 122 56\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"128 56 147 75 175 75\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"100 43 150 43 150 57 100 57 100 43\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"125.000000\" y=\"86.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"125.000000\" y=\"26.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_0</text>\n<polyline points=\"100 43 150 43 150 47 100 47 100 43\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"175 25 225 25\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"180 40 189 40 203 10 194 10 180 40 189 40 180 40\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"197.000000\" y=\"38.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=3.021948</text>\n<polyline points=\"75 125 103 125 122 144\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"128 144 147 125 175 125\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"75 175 103 175 122 156\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"128 156 147 175 175 175\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"100 143 150 143 150 157 100 157 100 143\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"125.000000\" y=\"186.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"125.000000\" y=\"126.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_1</text>\n<polyline points=\"100 143 150 143 150 147 100 147 100 143\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"175 125 225 125\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"180 140 189 140 203 110 194 110 180 140 189 140 180 140\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"197.000000\" y=\"138.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=0.494093</text>\n<polyline points=\"175 75 225 75\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"225 75 253 75 272 94\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"278 94 297 75 325 75\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"225 125 253 125 272 106\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"278 106 297 125 325 125\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"250 93 300 93 300 107 250 107 250 93\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"275.000000\" y=\"136.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"275.000000\" y=\"76.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_2</text>\n<polyline points=\"250 93 300 93 300 97 250 97 250 93\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"325 75 375 75\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"330 90 339 90 353 60 344 60 330 90 339 90 330 90\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"347.000000\" y=\"88.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=5.291995</text>\n<polyline points=\"75 225 175 225\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"175 175 203 175 222 194\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"228 194 247 175 275 175\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"175 225 203 225 222 206\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"228 206 247 225 275 225\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"200 193 250 193 250 207 200 207 200 193\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"225.000000\" y=\"236.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"225.000000\" y=\"176.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_3</text>\n<polyline points=\"200 193 250 193 250 197 200 197 200 193\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"275 175 325 175\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"280 190 289 190 303 160 294 160 280 190 289 190 280 190\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"297.000000\" y=\"188.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=6.094472</text>\n<polyline points=\"225 25 375 25\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"375 25 403 25 422 44\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"428 44 447 25 475 25\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"375 75 403 75 422 56\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"428 56 447 75 475 75\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"400 43 450 43 450 57 400 57 400 43\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"425.000000\" y=\"86.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"425.000000\" y=\"26.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_4</text>\n<polyline points=\"400 43 450 43 450 47 400 47 400 43\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"475 25 525 25\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"480 40 489 40 503 10 494 10 480 40 489 40 480 40\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"497.000000\" y=\"38.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=3.232733</text>\n<polyline points=\"325 125 353 125 372 144\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"378 144 397 125 425 125\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"325 175 353 175 372 156\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"378 156 397 175 425 175\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"350 143 400 143 400 157 350 157 350 143\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"375.000000\" y=\"186.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"375.000000\" y=\"126.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_5</text>\n<polyline points=\"350 143 400 143 400 147 350 147 350 143\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"425 125 475 125\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"430 140 439 140 453 110 444 110 430 140 439 140 430 140\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"447.000000\" y=\"138.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=4.438145</text>\n<polyline points=\"475 75 503 75 522 94\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"528 94 547 75 575 75\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"475 125 503 125 522 106\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"528 106 547 125 575 125\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"500 93 550 93 550 107 500 107 500 93\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"525.000000\" y=\"136.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"525.000000\" y=\"76.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_6</text>\n<polyline points=\"500 93 550 93 550 97 500 97 500 93\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"575 75 625 75\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"580 90 589 90 603 60 594 60 580 90 589 90 580 90\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"597.000000\" y=\"88.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=0.473583</text>\n<polyline points=\"275 225 425 225\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"425 175 453 175 472 194\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"478 194 497 175 525 175\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"425 225 453 225 472 206\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"478 206 497 225 525 225\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"450 193 500 193 500 207 450 207 450 193\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"475.000000\" y=\"236.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"475.000000\" y=\"176.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_7</text>\n<polyline points=\"450 193 500 193 500 197 450 197 450 193\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"525 175 575 175\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"530 190 539 190 553 160 544 160 530 190 539 190 530 190\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"547.000000\" y=\"188.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=1.289718</text>\n<polyline points=\"525 25 625 25\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"625 25 653 25 672 44\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"678 44 697 25 725 25\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"625 75 653 75 672 56\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"678 56 697 75 725 75\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"650 43 700 43 700 57 650 57 650 43\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"675.000000\" y=\"86.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"675.000000\" y=\"26.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_8</text>\n<polyline points=\"650 43 700 43 700 47 650 47 650 43\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"725 25 775 25\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"730 40 739 40 753 10 744 10 730 40 739 40 730 40\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"747.000000\" y=\"38.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=1.424995</text>\n<polyline points=\"575 125 603 125 622 144\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"628 144 647 125 675 125\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"575 175 603 175 622 156\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"628 156 647 175 675 175\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"600 143 650 143 650 157 600 157 600 143\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"625.000000\" y=\"186.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"625.000000\" y=\"126.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_9</text>\n<polyline points=\"600 143 650 143 650 147 600 147 600 143\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"675 125 725 125\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"680 140 689 140 703 110 694 110 680 140 689 140 680 140\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"697.000000\" y=\"138.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=2.780282</text>\n<polyline points=\"725 75 753 75 772 94\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"778 94 797 75 825 75\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"725 125 753 125 772 106\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"778 106 797 125 825 125\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"750 93 800 93 800 107 750 107 750 93\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"775.000000\" y=\"136.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"775.000000\" y=\"76.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_10</text>\n<polyline points=\"750 93 800 93 800 97 750 97 750 93\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"825 75 875 75\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"830 90 839 90 853 60 844 60 830 90 839 90 830 90\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"847.000000\" y=\"88.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=0.467264</text>\n<polyline points=\"525 225 675 225\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"675 175 703 175 722 194\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"728 194 747 175 775 175\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"675 225 703 225 722 206\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"728 206 747 225 775 225\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"700 193 750 193 750 207 700 207 700 193\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"725.000000\" y=\"236.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"725.000000\" y=\"176.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_11</text>\n<polyline points=\"700 193 750 193 750 197 700 197 700 193\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"775 175 825 175\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"780 190 789 190 803 160 794 160 780 190 789 190 780 190\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"797.000000\" y=\"188.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=4.883668</text>\n<polyline points=\"775 25 875 25\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"875 25 903 25 922 44\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"928 44 947 25 975 25\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"875 75 903 75 922 56\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"928 56 947 75 975 75\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"900 43 950 43 950 57 900 57 900 43\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"925.000000\" y=\"86.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"925.000000\" y=\"26.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_12</text>\n<polyline points=\"900 43 950 43 950 47 900 47 900 43\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"975 25 1025 25\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"980 40 989 40 1003 10 994 10 980 40 989 40 980 40\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"997.000000\" y=\"38.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=0.945475</text>\n<polyline points=\"825 125 853 125 872 144\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"878 144 897 125 925 125\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"825 175 853 175 872 156\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"878 156 897 175 925 175\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"850 143 900 143 900 157 850 157 850 143\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"875.000000\" y=\"186.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"875.000000\" y=\"126.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_13</text>\n<polyline points=\"850 143 900 143 900 147 850 147 850 143\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"925 125 975 125\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"930 140 939 140 953 110 944 110 930 140 939 140 930 140\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"947.000000\" y=\"138.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=3.37508</text>\n<polyline points=\"975 75 1003 75 1022 94\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1028 94 1047 75 1075 75\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"975 125 1003 125 1022 106\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1028 106 1047 125 1075 125\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1000 93 1050 93 1050 107 1000 107 1000 93\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"1025.000000\" y=\"136.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"1025.000000\" y=\"76.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_14</text>\n<polyline points=\"1000 93 1050 93 1050 97 1000 97 1000 93\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"1075 75 1125 75\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"1080 90 1089 90 1103 60 1094 60 1080 90 1089 90 1080 90\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"1097.000000\" y=\"88.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=1.600605</text>\n<polyline points=\"775 225 925 225\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"925 175 953 175 972 194\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"978 194 997 175 1025 175\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"925 225 953 225 972 206\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"978 206 997 225 1025 225\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"950 193 1000 193 1000 207 950 207 950 193\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"975.000000\" y=\"236.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"975.000000\" y=\"176.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_15</text>\n<polyline points=\"950 193 1000 193 1000 197 950 197 950 193\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"1025 175 1075 175\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"1030 190 1039 190 1053 160 1044 160 1030 190 1039 190 1030 190\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"1047.000000\" y=\"188.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=1.029238</text>\n<polyline points=\"1025 25 1125 25\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"1125 25 1153 25 1172 44\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1178 44 1197 25 1225 25\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1125 75 1153 75 1172 56\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1178 56 1197 75 1225 75\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1150 43 1200 43 1200 57 1150 57 1150 43\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"1175.000000\" y=\"86.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"1175.000000\" y=\"26.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_16</text>\n<polyline points=\"1150 43 1200 43 1200 47 1150 47 1150 43\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"1225 25 1275 25\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"1230 40 1239 40 1253 10 1244 10 1230 40 1239 40 1230 40\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"1247.000000\" y=\"38.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=3.533838</text>\n<polyline points=\"1075 125 1103 125 1122 144\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1128 144 1147 125 1175 125\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1075 175 1103 175 1122 156\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1128 156 1147 175 1175 175\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1100 143 1150 143 1150 157 1100 157 1100 143\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"1125.000000\" y=\"186.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"1125.000000\" y=\"126.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_17</text>\n<polyline points=\"1100 143 1150 143 1150 147 1100 147 1100 143\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"1175 125 1225 125\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"1180 140 1189 140 1203 110 1194 110 1180 140 1189 140 1180 140\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"1197.000000\" y=\"138.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=1.737775</text>\n<polyline points=\"1225 75 1253 75 1272 94\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1278 94 1297 75 1325 75\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1225 125 1253 125 1272 106\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1278 106 1297 125 1325 125\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1250 93 1300 93 1300 107 1250 107 1250 93\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"1275.000000\" y=\"136.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"1275.000000\" y=\"76.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_18</text>\n<polyline points=\"1250 93 1300 93 1300 97 1250 97 1250 93\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"1325 75 1375 75\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"1330 90 1339 90 1353 60 1344 60 1330 90 1339 90 1330 90\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"1347.000000\" y=\"88.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=0.632872</text>\n<polyline points=\"1025 225 1175 225\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"1175 175 1203 175 1222 194\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1228 194 1247 175 1275 175\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1175 225 1203 225 1222 206\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1228 206 1247 225 1275 225\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"round\" />\n<polyline points=\"1200 193 1250 193 1250 207 1200 207 1200 193\" fill=\"black\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"1225.000000\" y=\"236.000000\" font-size=\"7.000000\" text-anchor=\"middle\"></text>\n<text x=\"1225.000000\" y=\"176.000000\" font-size=\"7.000000\" text-anchor=\"middle\">theta=theta_19</text>\n<polyline points=\"1200 193 1250 193 1250 197 1200 197 1200 193\" fill=\"lightgray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<polyline points=\"1275 175 1325 175\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"1280 190 1289 190 1303 160 1294 160 1280 190 1289 190 1280 190\" fill=\"gray\" stroke=\"black\" stroke-width=\"1.000000\" stroke_linejoin=\"miter\" />\n<text x=\"1297.000000\" y=\"188.000000\" font-size=\"7.000000\" text-anchor=\"start\">Φ=5.036813</text>\n<polyline points=\"1275 25 1375 25\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"1325 125 1375 125\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"1325 175 1375 175\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<polyline points=\"1275 225 1375 225\" fill=\"transparent\"stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\" />\n<path d=\"M 1375 25 L 1390 25\" fill=\"none\" stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\"/>\n<text x=\"1400.000000\" y=\"25.000000\" font-size=\"12.000000\" text-anchor=\"end\">0</text>\n<path d=\"M 1375 75 L 1390 75\" fill=\"none\" stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\"/>\n<text x=\"1400.000000\" y=\"75.000000\" font-size=\"12.000000\" text-anchor=\"end\">1</text>\n<path d=\"M 1375 125 L 1390 125\" fill=\"none\" stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\"/>\n<text x=\"1400.000000\" y=\"125.000000\" font-size=\"12.000000\" text-anchor=\"end\">2</text>\n<path d=\"M 1375 175 L 1390 175\" fill=\"none\" stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\"/>\n<text x=\"1400.000000\" y=\"175.000000\" font-size=\"12.000000\" text-anchor=\"end\">3</text>\n<path d=\"M 1375 225 L 1390 225\" fill=\"none\" stroke=\"darkred\" stroke-width=\"3.000000\" stroke-linejoin=\"miter\"/>\n<text x=\"1400.000000\" y=\"225.000000\" font-size=\"12.000000\" text-anchor=\"end\">4</text></svg>"
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"bs = pcvl.Circuit.generic_interferometer(N,\n",
" lambda idx : phys.BS(theta=pcvl.P(\"theta_%d\"%idx))//(0, phys.PS(phi=np.pi*2*random.random())),\n",
" shape=\"rectangle\",\n",
" depth = 2*N,\n",
" phase_shifter_fun_gen=lambda idx: phys.PS(phi=np.pi*2*random.random()))\n",
"pcvl.pdisplay(bs, recursive = True)\n"
]
},
{
"cell_type": "markdown",
"id": "6059f2d3",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"### Minimizing infidelity"
]
},
{
"cell_type": "markdown",
"id": "8bb792d0",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"We first define a \"loss\" function which computes the infidelity between the target unitary and the implemented unitary. The fidelity evaluates the performance of the multiport interferometers via $F=\\frac{\\lvert \\mathrm{Tr} ({U_{target}}^{\\dagger} U_{interf}) \\rvert²}{N\\mathrm{Tr}({U_{target}}^{\\dagger} U_{target})}$ where $U_{target}$ is the target matrix and $U_{interf}$ the actual unitary matrix realized by the interferometer and where N is the size of the matrices."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "2faecde6",
"metadata": {
"pycharm": {
"is_executing": true,
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"def infidelity(c, U_target, params, params_value):\n",
" \"\"\"\n",
" inputs: \n",
" *c* type: circuit from pcvl.Circuit()\n",
" role: generic interferometer to optimize \n",
" *U_target* type: matrix from pcvl.Matrix.random_unitary()\n",
" role: unitary matrix randomly chosen\n",
" *params* type: parameters from c.get_parameters()\n",
" role: get the parameters associated to the circuit c \n",
" *params_value* type: int\n",
" role: value of the parameters\n",
"\n",
" outputs: \n",
" *infidelity value* type: int\n",
" role: infidelity value between U and U0\n",
" \"\"\"\n",
" for idx, p in enumerate(params_value):\n",
" params[idx].set_value(p) #give a value to each params\n",
" U = c.compute_unitary(use_symbolic=False) #\n",
" U_dag = np.transpose(np.conjugate(U))\n",
" f = abs(np.trace(U_dag @ U_target)) ** 2 / (c.m * np.trace(U_dag @ U))\n",
" return 1 - abs(f)"
]
},
{
"cell_type": "markdown",
"id": "81adaa92",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"### Basinhopping algorithm"
]
},
{
"cell_type": "markdown",
"id": "fa4c270f",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"Then, we want to optimise over the angles of the phase-shifters $\\phi$ comprising $U_{interf}$ by minimizing the infidelity $1 - F$ [1]. To do so, we use a numerical optimisation algorithm based on the basinhopping algorithm to explore the space of phases and minimize $1-F$ for a specific value of manufacturing error $\\alpha$. We then vary $\\alpha$ to see how robust to manufacturing errors our scheme is."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "49aecaa7",
"metadata": {
"pycharm": {
"name": "#%%\n",
"is_executing": true
}
},
"outputs": [],
"source": [
"def calculate_angle(index_unitary,U_target,alpha):\n",
" \"\"\"\n",
" inputs:\n",
" *index_unitary* type: int\n",
" role: index on the number of unitary matrices \n",
" *U_target* type: matrix from pcvl.Matrix.random_unitary()\n",
" role: unitary matrix randomly chosen\n",
" *alpha* type: int\n",
" role: angle error due to manufacturing \n",
"\n",
" outputs:\n",
" *infidelity min* type: int\n",
" role: min value of infidelity\n",
" \"\"\"\n",
" start = time.time()\n",
" c = pcvl.Circuit.generic_interferometer(N,\n",
" lambda idx: (phys.BS(theta=(45+alpha)/180*np.pi)\n",
" // (0, phys.PS(phi=pcvl.P(\"phi_m%d\" % idx)))),\n",
" depth=2*N,\n",
" phase_shifter_fun_gen=lambda idx: phys.PS(phi=pcvl.P(\"phi_r%d\" % idx)))\n",
" params = c.get_parameters() #We get the parameters of the circuit we have just created\n",
" \n",
" infidelities = []\n",
"\n",
" for _ in range(n_try): #We do a number of tests\n",
" init_params = np.random.randn(len(params)) \n",
" res = basinhopping(lambda x: infidelity(c, U_target, params, x), init_params, stepsize=0.1, niter=n_iter) #the algorithm is looking for a global minimum of infidelity\n",
" infidelities.append(res.fun) #fun is the value of the function at the solution\n",
" return min(infidelities)"
]
},
{
"cell_type": "markdown",
"id": "a2b378aa",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"Based on the definition of the loss function, we now explore the space of phases to obtain a correct implementation of the target unitary with our BS-based scheme. We then reiterate the following function for various error parameter $\\alpha$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f12a5609",
"metadata": {
"pycharm": {
"is_executing": true,
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"n_angles = (abs(angle_min)+abs(angle_max))/angle_step\n",
"\n",
"def discovery_unitary(index_unitary):\n",
" \"\"\"\n",
" inputs:\n",
" *param index_unitary* type: int\n",
" role: index on the number of unitary matrices\n",
" outputs:\n",
" type: list of int\n",
" role: infidelities for all alpha angles\n",
" \"\"\"\n",
" np.random.seed(int(10000*time.time()+os.getpid()) & 0xffffffff) #not necessary here, but when we run our program on several cores, this line allows us to generate different unit matrices at the same time.\n",
" U_target = pcvl.Matrix.random_unitary(N) #generates the unit random matrix\n",
" l_infidelities = []\n",
" M=angle_min + math.ceil((angle_max-angle_min)/angle_step)*angle_step + 1\n",
"\n",
" for alpha in tqdm_notebook(range(angle_min,M,angle_step), desc = \"Angle progress bar unitary {}\".format(index_unitary+1), leave=False):\n",
" l_infidelities.append(calculate_angle(index_unitary, U_target, alpha)) \n",
" return l_infidelities\n",
" \n",
"try:\n",
" os.remove(\"%s-%d.log\" % (logfilebs,N)) #To avoid having data on our file that we don't want, we delete the data from previous runs\n",
"except OSError:\n",
" pass\n",
"\n",
"for index_unitary in tqdm_notebook(range(n_unitary), desc = \"Unitary progress bar\"):\n",
" result = discovery_unitary(index_unitary)\n",
" e = datetime.datetime.now()\n",
" with open(\"%s-%d.log\" % (logfilebs,N), \"a\") as f: #we create our file\n",
" f.write(e.strftime(\"%Y-%m-%d %H:%M:%S\")+\"\\t\") #we add the date and the hour of calculation \n",
" f.write(\"\\t\".join([\"%g\"%inf for inf in result])) #we add the list of infidelities returned by discovery_unitary()\n",
" f.write(\"\\n\")"
]
},
{
"cell_type": "markdown",
"id": "c42f7fe3",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"## Result of the BS-based scheme"
]
},
{
"cell_type": "markdown",
"id": "a3f506da",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"In this section, we treat the previously collected infidelities as a function of $\\alpha$."
]
},
{
"cell_type": "markdown",
"id": "9c4f8e3f",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"The collected data looks like that the output of the following cell. Each lign starts with the date and hour of collection and then, we have values for the minimal infidelity with several values of $\\alpha$. We have a lign per unitary matrix $U_{target}$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "1fd3af77",
"metadata": {
"pycharm": {
"is_executing": true,
"name": "#%%\n"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"2022-06-20 10:53:31\t0.0358288\t2.10034e-11\t6.11298e-11\t0.0214913\t0.181749\t0.575823\n",
"\n",
"2022-06-20 10:54:00\t5.058e-11\t6.30681e-11\t5.92989e-11\t5.84721e-12\t0.00633146\t0.24023\n",
"\n",
"2022-06-20 10:54:41\t5.27425e-11\t1.92772e-11\t1.84568e-11\t0.0060772\t0.0685131\t0.344481\n",
"\n"
]
}
],
"source": [
"with open(\"bsbasednotebook-opt-5.log\", \"r\") as filin:\n",
"... for ligne in filin:\n",
"... print(ligne)"
]
},
{
"cell_type": "markdown",
"id": "7a16480e",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"### Graph of the infidelity as a function of $\\alpha$ in degree : BS-based"
]
},
{
"cell_type": "markdown",
"id": "25232b73",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"We plot the infidelity as a function of $\\alpha$. To do this, we compute the average of the infidelities over all $U_{target}$ tried."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b8168095",
"metadata": {
"pycharm": {
"is_executing": true,
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"L=[]\n",
"with open(\"bsbasednotebook-opt-5.log\", \"r\") as flog:\n",
" for l in flog:\n",
" L.append([float (f) for f in l.strip().split(\"\\t\")[1:]])\n",
"\n",
"A=np.asarray(L)\n",
"print(A)\n",
"print(A.shape)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4ce9d2d6",
"metadata": {
"pycharm": {
"is_executing": true,
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"Lmin=[]\n",
"Lmax=[]\n",
"Laverage=[]\n",
"\n",
"for i in range(A.shape[1]):\n",
" Langle = A[:,i]\n",
" Laverage.append(np.average(Langle))\n",
" Langle.sort()\n",
" Lmin.append(np.average(Langle[0:10]))\n",
" Lmax.append(np.average(Langle[-10:]))\n",
"\n",
"print(Laverage)\n",
"print(Lmin)\n",
"print(Lmax)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "6d77b08a",
"metadata": {
"pycharm": {
"is_executing": true,
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, 'infidelity 1-F')"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEgCAYAAABM0P/cAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8bklEQVR4nO3dd5gUZbb48e/pnpwzORhAzKKIiglUFFwVs5gREHdX13VNd+/d367evWHvml11VZIYUcyoKCYQMIKKARFFEQGByTn1dJ/fH9UjvePMMMOE6nA+z1NPT79dXX2mYOp0vVFUFWOMMaYzPG4HYIwxJvJY8jDGGNNpljyMMcZ0miUPY4wxnWbJwxhjTKdZ8jDGGNNpljyM6WYislREorIPvIjMExEVkaFux2LcZcnDRJXgha3l1iAiP4jIwyKydxvvO1REHheRjcH9K0XkOxF5SURuFJHU3v5djAlncW4HYEwP+c+QnzOB0cAlwFkicpSqrm5+UUQuAh4GBHgbeB6oA4YARwGnAM8B63slcmMigCUPE5VU9eaWZSJyD3AVcA0wJViWAtwHKHCiqr7VyvvGAMU9F60xkceqrUwseT34mB9Sth+QAXzZWuIAUNX3VLW8sx8mIoki8t8isiFYFfadiNwkIgmt7Hu6iDwmIt+ISE1w+1hErhaRX/ydikgfEblNRNYF9y0P/jxPRHZvZf+TRGSRiBSHxHKriGS1EfsJIrI8eOxSEXlBREZ09hyY6GV3HiaWnBB8XBVSVhJ87C8iqapa042ftwA4FHgG8AGTgJuBUSJymv7rxHL/BwSAD4EtOFVtxwF3B49xcfOOwbuld4E9gDeAl3Cq3IYEP+MZ4PuQ/W8Kfm4p8DJQCBwAXA+cLCJHqGplyP5nA08BjcHHrTjVd+8Dn3f1pJgooaq22RY1G071k+JcLJu3O4DlOBfnl4D0kP0F+Cj4ntXAlcBIIKELMSwNHu8bIDukPAnnAqzAxS3es0crx/HgtMUocFhI+anBsjtbeU9Ci99vXHDf94CsFvtOaXkcIA0nofqAUS32vzPk/A51+9/aNnc31wOwzbbu3EIubq1ta4ALWnnPYGBJi30bce4C/g3I6GQMzcnj4lZeGxt8bUkHj3VwcP+/hJQ1J4//7cD7nw/uu28br38KFIY8vzC4/8Ot7JsJlFvysE1VrdrKRCdVleafg91s98WpGnpcRPZV1T+F7PsjMC7YjXc8MAqnd1bz9lsRGauqG4LHuwbIavGRL2hID66gd1oJbQXgx7m7+ZmI5AI3ACcDuwMtuwYPaHHcLcAfReRgYBFONdZqVfW3eN8ROHcR54jIOa3EkwDki0iuqpbgJKtWY1fVChFZDRzbynFMjLHkYaKeOu0YH4nImcBm4EYReUBVN7XYby2wtvl5sIF4Ls4F+E7g9OBL1+C0L4T6AafaK9T2VmJpEpFioCDkc7KAlcBuOFVoj+C0TzThJKnfA4khx6gUkcNxuiOfBpwUfKlYRP4J/Leq+oJluTh/5ze1jKWF5uqqzLZiD9q2k+OYGGHJw8QMVS0XkXU4364PBjbtZP+vReRinPEdx4WUD+3gR/YBfgwtEJE4IA+oDCmejpM4/lNbdDEWkSNwkkfL2DYD00REgH2C8V0J/AWnreTPwV0rAI+q5nQw5oqQ2FvTt4PHMVHOuuqaWJMdfOzo//2q4KO0u1frWqveOQrw4rQ1NNsz+PhsB4/xM3WsUdV7cKrcYMcdEsAHQLaI7NuhiOGTtj5XRDKBgzp4HBPlLHmYmCEip+N8w/fh9D5CRHYLjqXIbGV/AZrbRpbtwkf+WUSakxUikgT8Lfj0oZD9fgg+jm3x+SOBf28lrn1FpLU7g+ay2pCyO4OPs0SkfyvHSg1WgTV7ESgDLhCRUS12v5kd1Vomxlm1lYlKInJzyNNUnKqdicHn/6GqzXX6mThjKW4VkXeBL3HuNgpwqoJ2xxkXcd0uhLEWWCMioeM89gBeAR4N2e8RnMbyu0RkHPAtMIwd06Kc1+K444Pxvo/THbgQGBg8fgC4tXlHVX1LRP6Ik7S+FZFFwAacNo4hOHcYK4AJwf2rRWQGzviO5SISOs5jP5wkeswunAsTbdzu7mWbbd250XoX3SacC+CLwPgW+yfiVPP8E6fKZjvOhb4C+Bj4byC/kzEsDX5uYvD9G4AGnIF7NwGJrbxnH2AhTiKoCX72dGBo8FjzQvbdG2fsyiqgKHjsH3AGB45pI6ajcAYt/oTTDbkIp4H/DlqM5wjuPx4nqdTi3Im8CIwA5mFddW1TRVSjcuZoY4wxPcjaPIwxxnSaJQ9jjDGdZsnDGGNMp1nyMMYY02kx01U3Ly9Phw4d6nYYxhgTVhqbAlTW+qmsUmprobHeg/q88PP0cB8Xq2p+y/fFTPIYOnQoq1at2vmOxhgTpYqrG/hgXQWL3q7ngw+F775KoGZzBv6qZGcHT4CE/Crydq9mvwP9HHuklz9dPHBja8eKmeRhjDGxpK7Rz+qNFSx6p5Zl7wX4+vM4yjam4yvOp3m2nbisGjKGVjB8vyLGHCH8amwyo/bMJDN5x0QCf7q49eNb8jDGmAjnDyjfbK/ijQ+reGtZE599Ese29ck0bs9Em5w5MT3JDSQPqGD44eWMOlSZMDaBYw7IpF9mH5yZeDrHkocxxkQQVeWninqWfV7Ba0sbWfWRsHFdEnVbMgnUZQAgcX4S+lYw5MitHHBwgBOO8TL+sAyG9cnH69mVOT5/yZKHMcaEsYpaHyu/K2fR0nree1/55ssEqjal01TePDu+Ep9XRf6+xYw4oIljxniYeEwKBw3JJDmhozPxd54lD2OMCRMNTX6+3FzF4hU1LF3hZ81ncRRvSKWxMBfUGVnhTa8jbVAle5xQwmGHCSePS+SIvTLJTcvo1VgteRhjjAsCAeX74hqWfFzF60t9rP7Ew5Zvkqnflok2ZgEgiT6S+lcw/MQtHHyIcuLYeI47OIPBOQW71E7RnSx5GGNMLyisque9rypZ9HYDH30E369NoHZLJv7qfs4OngAJfSoZcOh29jvIz7ijvUwYk8aIfjkkxIXfeG5LHsYY081qGpr4ZEMFr75Tx4r3A3z9eTxlG9NpKv156XricqrJ2rOMvfbzceQY4eRjkjlkj0zSk7LcC7wTIjJ5iMjuOCu8Zarq2W7HY4yJXU3+AF9vq+L192p4e3kTX3zqZfv3KTRszwJ/LgCelAZSB1Yw9OgyRo9WTjo2kaP3y6RvZj93g++CXk8eIjIXZ4W0QlXdL6R8As6Kbl5gtqr+X1vHUNXvgWnBFdqMMaZXqCqby+pYsrqS15c28vEq4cd1SdT/lEmg3hlYJ/FNJPatYPdjt3LQwX5OODae40els2dBPp5u6iYbDty485gH3Iuz9CYAIuIF7sNZvWwzsFJEFuIkkr+1eP9UVS3snVCNMbGsrKaRD74p59UlDbz/obJ+TQJVmzLwVwa7yYoSn19JnwOK2ftAp5vshKNTOXBwFknxue4G38N6PXmo6jIRGdqieDSwPnhHgYg8CUxS1b/h3KXskuBazDMABg8evKuHMcbEiB9LannstTLeWurnq8/jKNmQ5kznEZwk0JtZS/rgSobtW8zhhwkTxyZy+PAsslMzd3Lk6BMubR4DgE0hzzcDh7W1s4jkAv8DjBSRfw8mmV9Q1ZnATIBRo0bZervGmH+hqny5uYpZz5bz8kIPmz7Npql8AACepEaS+lew56EVHHKIcuK4BMYemMHA7F2bziPahEvy6BRVLQF+7XYcxpjI0+QP8O43ZcycX8Obr8ZR/FUugdrB4PWTtlsph5xRwmknx3HiEamM6JtLnDf8usmGg3BJHluAQSHPBwbLjDGmy+p9fhatKmH2E3W893YSld/mor5cJMFH1l4lHDO+kekXJDP+wBwS47xuhxsRwiV5rASGichuOEljMnCBuyEZYyJZeW0jC5aW8MhTjXy6PJXajXkQ8OBNq6ffqO2c9Cs/M85NZ/QefaKqF1RvcaOr7nxgLJAnIpuBm1R1johcBSzG6WE1V1XX9HZsxpjItrmsjodfKeGpZwJ882EGDduccRRxOdXsefwWTj8dpp+RxfC+/a3doovc6G11fhvli4BFvRyOMSaCqSprt1Yx65kKFi6EHz/JpqlsIACJ/csYedaPnHeWl0tOzqFf5qCdHM10RrhUWxljTIf4A8p735Qx68lq3ng1jsI1OQRqBoEnQOpuJYw5vYRLJydw1lF5ZKZkux1u1LLkYYwJe/U+P4s/LWH2E/WseDOBim/z0MYcp8F7eAlHndDA9AuSOfGgXJLi890ONyZY8jDGhKWKWh9PLy/ikfk+PlmeSs0PToO3J7WePiO3c9LJfmacl8Zhw/p02+p4puMseRhjwsbWijoeXlTKk08H+PqDdBq29gcgLruaPcb9xGmTlMvPymJEP2vwdpslD2OMa1SVdduqmfVsOQtfhB8+yaap1BnhndivnAPO+JHJZ3u5eGIOA7MHuhytCWXJwxjTq/wBZeV3ZTwwv5rXF3kp/DIXf7DBO2VIKYefUsolk+M559hcslKy3A7XtMGShzGmx9X7/Lz5WSmz59ey7I1EytflBhu8m8gcVsyRxzcw/fxkJhySS1J8ntvhmg6w5GGM6RGV9T6eXV7Mw/N9rFqWTM2GXAjk40lpoODAIsaf3MTl56Vy5F7W4B2JLHkYY7rNtop6HllcwpML/Hz1QToNW4IjvLNq2G3sT5x6qnL52ZnsO6CfNXhHOEsexpgu+WZbFbOeK+fFF+GHj7PwlTgN3gl9y9l/0o+cc5aHS0/OYXCuNXhHE0sexphOCQSUld+XM/OpKl57xcv2L3LxVw8CcRq8D51YysWTEzh3bC45qVluh2t6SMwkj6aArQVlzK5qaPLz1melzHqylmWvJ1C2Lg9tyEbim8gYVsKY47Yz9fwkfjUqj+QEa/COBTGTPL7Z4HM7BGMiSlW9j+ffK+GhJxpY9U4K1RtywJ+PJ7mB/P2LOH6CjxmT0zhqRL4tmBSDYiZ51FXFsbG4liF5KW6HYkzYKqys59HXS5i/wM+a99Oo39wHEOIyaxl69FZOOVW5/JwM9h9oDd6xLmaSB34Pd80v4c7fWfIwpqXbF/zEzIf8bFiVha842ODdp4J9T93M2WfCZafkMiTPGrzNDqIaG20BIqO0/xEvsvldmxPHmFAvrSxi0lHZqM9D8uAyDjyqhovOi2PyuFxy0xLdDs+4TEQ+VtVRLctj5s7Dm+xj26f5rPyunNF72hz/xjS7/f46tDGf6bds4O7fDyIlIdftkEwEiJlWrswsJVCfwO1zK90OxZiwUV7byIevpxOfV8VfL+9HSkLMfJ80XRQzyaNPnhdPSgOLX0zC5w+4HY4xYWHmS0XUb8lm/3Gl9MtKcjscE0FiJnmkJHgZcEgRFevyeGVlsdvhGBMWZs8NgChXXZ7gdigmwsRM8gA4/0IFv5e7Z9e5HYoxrlu3tYrv388jbY9izh9b4HY4JsLEVPK4+tw84nOr+eD1dCrqbNCgiW23zC3DX5XM8afVkhTvdTscE2FiKnkMyE5mxFGl1G/K4ZHFhW6HY4xr/AHlxafj8ST6uGFGutvhmAgUU8kDYMZUpzfJzIf8LkdijHsWry6hdE0BfUYWMma4dV03nRdzyePSE/NJHlzK1yty2FRa63Y4xrji9gdq0CYv51/ot0GzZpfEXPJIT4rn8BOraCpN456nStwOx5heV1nv4/3X0onLqeYPF+a7HY6JUDGXPACuuTwZvH7mPyHEyvQsxjSb80oRdZty2PfYUgZmJ7sdjolQMZk8Jh6SR9aIIn76JJ9Pf7AR5ya2zJzjB5QrL493OxQTwWIyecR7PUyY1EigNpHbHyp3Oxxjes13hTV8+24uqbuXcOHxVmVldl1MJg+A66Zm4Elu5NXnE2my6UpMjLh1Xgn+yhTGnlJj81iZLonZ5HHI7pn0G1lE+df5vPqJNZyb6BcIKM89FYck+Lh+Rprb4ZgIF7PJQ0SYfIEfbfJy12zrsmui31tflFD8RQEFBxZxzN45bodjIlzMJg+AqyfnEZddzfuL06iqt+lKTHS77cEa1BfHuRc04fHY2A7TNRGZPERkbxF5QESeEZHf7OpxBuemsNeRZdRtzOWxN4u6M0Rjwkp1QxMrFqURl13DdRfluR2OiQK9njxEZK6IFIrIly3KJ4jIOhFZLyJ/bO8YqrpWVX8NnAsc2ZV4pk91TsGDc5u6chhjwtq8xYXUbsxlr6NLGJKX4nY4Jgq4cecxD5gQWiAiXuA+YCKwD3C+iOwjIvuLyMsttoLge04DXgEWdSWYyyYUkDSwlK+WZ/NTuU3VbqLTA7OdL0e/mW49rEz36PXkoarLgNIWxaOB9ar6vao2Ak8Ck1T1C1U9pcVWGDzOQlWdCFzY1meJyAwRWSUiq4qKWq+WykyOZ/T4KnzF6dz7tC0SZaLPxuJa1i3PJWVoCZeeaOt2mO4RLm0eA4BNIc83B8taJSJjReQfIvIg7dx5qOpMVR2lqqPy89seEHX1dGe6kicet+lKTPS57dFimspTOfrkatIS7c7DdI+I/J+kqkuBpd11vFNH55E5vIhNq/L4fFMlBw7O7K5DG+OqQEB5en4cktDE9Vekuh2OiSLhcuexBRgU8nxgsKxXJMR5GH9aPYGaJO6w6UpMFFm2tpTCz/LJ26+Q4/bLdTscE0XCJXmsBIaJyG4ikgBMBhb2ZgDXTcvAk9TIy88l2HQlJmrc+mA12hjPWZNtbIfpXm501Z0PvA/sJSKbRWSaqjYBVwGLgbXAAlVd05txHbZnFn0OLKLsqwLe/Kxle74xkae2sYl3XknFm1nL9ZfaXYfpXm70tjpfVfuparyqDlTVOcHyRao6XFX3UNX/6e24RISzJzvTldw5p6a3P96YbvfYm0XUfJ/LsCNL2KPA2jtM9wqXaquwcO1FucRl1bDi1VRqGmzQoIls9832AcIV0+3P3HS/Nv9XiUh/EYnI3li7amheKnuOKaV2Qy5PLCl0Oxxjdtmm0lq+eieH5EGlTJ1oYztM92vvK8km4ODmJ+J4RESG9HxY7pk2xQMI98+xOw8Tue58vISm0jTGTKwiI8lWDDTdr73k0bJrhge4CIjqlrfpv8onsX8Za97JorCy3u1wjOk0VeXJxz1InJ9rbWyH6SFWGdpCVkoChxxfSWNRBvc9a9OVmMjz3jdlbFtdQM5+hZx0UFR/1zMusuTRiqsvTwRPgEcedTsSYzrv1plVaEM8p5/jw2tjO0wP2VmD+BgRaZ783wMocKSI9G25o6p2aXbbcHL6EflkDCtm08o81mypZN8BGW6HZEyH1Pv8vLUwBW96HTdclu12OCaK7Sx53NFK2d2tlCng7Xo44SExzstxp9bxwm0F3D6vkLl/suRhIsP8JYVUf9eXYeM3s1e/QTt/gzG7qL3ksVuvRRGGrpuezov3+Hjp2Xj8/652+28iwj2zfKDC9GlWI216VpvJQ1U39mYg4ebI4dkUHLCVwtUFLPmyhBMOsKU7TXjbWl7Pl0uzSRpQxoxT2l6CwJju0OmvJ8HxHnNFZHBPBBQuRISzzmtCfXHcMcumKzHh766nivAVp3PYSZVkpSS4HY6Jcrtyb+sBLgWi/qv4Hy7KwZtRy7JXU6hr9LsdjjFtUlUef8QDXj/XXpHsdjgmBuxqxWhMNADs2SeN3Q8voeb7PJ5aZtOVmPD10fpytn6ST/Y+RZx8SNR/rzNhwFrVduKyKR5Q4b45jW6HYkybbp1dSaA+gVPPbiDOa3/Wpud1+n+ZqvqBy4AN3R9O+LnitHwS+5bz+ZIsiqoa3A7HmF+o9/l5/YVkvGn13DjNxnaY3rFLX1FU9WFVLevuYMJRTmoCI4+voHF7Jve/UOR2OMb8wjMriqn6No8hhxXZgFbTa7p8fysie4vIX7ojmHB11bQkkACPPKJuh2LML/xjZj2oh8umxERTpAkT3VE5ug9wUzccJ2yddVQeaXuU8MOHeXyzrcrtcIz5WWFVPZ+9nU1iv3J+e4at22F6T3uLQQ3uyAZE/WikpHgv406pxV+VzO0Px0RtnYkQ9y4oprEwg1HjK8hJtbEdpve0Nz3JDzhzVu2MdHC/iHbt5Wm8/E8fLzwdx/03KB6brsS4TFV5+GHA6+eaK5LcDsfEmPaSRxXwNjB7J8c4Crix2yIKU8eMyCF//60UfZHPsrWljN3X1kkw7lr9QyWbV+WTNaKY00ZHfQWACTPtJY+PgExVfaW9A4hITHzl8XiESWc3MevjeO6YXc3YOy15GHf9fXYFgbpMJp5ZT0Kcje0wvau9/3HLgGEdOEZRcN+od/2UbLzpdSx9JZl6n01XYtzT2BTgtecT8aTWc8PULLfDMTGozeShqv+lqjtdEEBVl6nquO4NKzwN75vO0MOKqVqfxzMrbMyHcc/z7xdRsS6fQYcWc9AQG9thep/d63bSpZcIqId7Z9toc+Oeu2c2QMDDlEudGaCN6W2WPDrp16fnk9CngtVvZ1JSbQnE9L6S6gY+eTOThD4VXHmWNZQbd1jy6KT89EQOGFtBw7YsZr5U7HY4Jgbd91wRDdsyGXl8OfnpiW6HY2KUJY9d8Ntp8SDKQw8H3A7FxKB58wBPgGtmxERHRxOmLHnsgsnHFpC2ewnfv5/Lt9ur3Q7HxJDPf6xk40d5ZAwv4owjbd0O454OJQ8RsUENIZITvBxzcg3+yhTueqzU7XBMDLl1bjmBmiROOr2exDiv2+GYGNbRO4+fRGSBiEwUEbtbAf4wPRWJb+LZBV5Uo352FhMGfP4ArzyTgCe5gRunZ7odjolxHU0EVwAFwMvAJhH5XxEZ3nNhhb9x++WSu28RhZ8V8O43Nlmi6Xkvf1RM2df5DBhVxCG7W/Iw7upQ8lDVeao6FmfE+RzgAmCtiLwrItNEJK0HYwxLXo9w6tk+tCGe22fbNO2m5935YB34vVx8idrYDuO6TlVBqer3qvoXVR0KjAf8wExgm4jME5GDeyDGXxCRsSKyXEQeEJGxvfGZrbl+SjbetHreesmmKzE9q6ymkY/eyCA+v5LfnWNjO4z7Ot1+ISIpIjIF+AvOjLpfAXcCewMrReSGnbx/rogUisiXLconiMg6EVkvIn/cSRgKVANJwObO/g7dZZ8B6Qw+1FkC9MUPbLoS03MeXFhEw0/ZHDiunL6Z1kXXuK/DyUNEjhGRh4BtwN3AOuBwVd1fVf+sqocB/w7s7MI/D5jQ4the4D5gIs7KhOeLyD4isr+IvNxiKwCWq+pE4N+A/+zo79ATLroICHi4x6YrMT1ozkMBkAC/u8IWfDLhoaNddb8DlgB7AlcD/VT1ClX9qMWubwHZ7R1LVZcBLfu3jgbWB6vFGoEngUmq+oWqntJiK1TV5tF5ZYCrQ2yvOjuP+PxKVr2VQVlNo5uhmCi19qcqNryfT/qexZx7tFVZmfDQ0TuPZ4B9VPXoYON5bWs7qerHqrorXXkHAJtCnm8OlrVKRM4UkQeBR4F729lvhoisEpFVRUU9U61UkJHEfseW07Alm9mLrOrKdL9b55bhr07ihEl1JMXb2A4THjp6oV8DtDqRk4jkiMgl3RfSzqnqc8E7n/NUdWk7+81U1VGqOio/v+e+sf12ehygzJ1n05WY7tXkD/Di0wl4khq54XKbet2Ej44mj4eAPdp4bbfg612xBQhdO2RgsCwinD+2gJShJax/L4cfimvcDsdEkVc/LaHsq3z6jSzi8GFZbodjzM86mjza61SeC1R2MY6VwDAR2U1EEoDJwMIuHrPXpCbGcdTEGprKU7nz8RK3wzFR5I4HatEmL+dfFLCxHSastLmGuYhMAiaFFP1ZRFpW6icBR+Nc/DtEROYDY4E8EdkM3KSqc0TkKmAx4AXmquqajh4zHPxheipvzPLz9JNe7rraBnGZrquo9fHB4nTic6u45nybBNGElzaTB850JPuHPN8D6Ntin0bgdeC/O/qBqnp+G+WLgEUdPU64GX9QLtn7bGP76nw++LacI4a32+nMmJ2a/Uoh9ZsHcNCZPzIge7Db4RjzL9pbw3yWqh6qqocC7wBnNT8P2Y5U1WmquqH3Qg5PXo9wypk+AvUJ3DG3q7V4xsDMuQEQ5aoZ8W6HYswvdHRuq3Gqurang4l0116WiSe1njcWJtHYZD2vzK77dns1372XS9ruJVwwrsDtcIz5hfbaPH4LPK2qRcGf26Oqen/3hhZ5DhiUwcBDtvDju/146aMizhrTx+2QTIS69aFS/JWDGTe1iOQEa+8w4ae9No97gVVAEe0MxAtSIOaTh4hw4UXK35Z5+cfses4a43ZEJhL5A8rzC+KRRB/Xz0h3OxxjWtVem4enefqR4M/tbTbsNeiqc/KIz6viozfSqajzuR2OiUBvflZCyZf59DmwkKNHWMcLE55sVcBu1j8rmb2PLqN+cw4PvVbodjgmAt32YA3qi+O8C2xshwlf7bV57NOZA6nqV10PJzr8emocv30eZj/k55oz3I7GRJKqeh/vvppGXHY1116U63Y4xrSpvTaPL3HaMnZGgvtZ1VXQxeMLuG5ICetW5PBjSS2Dc1PcDslEiHmvFVL34wD2n/Qjg3NtbIcJX+0lj3G9FkWUSUuMY8xJ1bw1M5e75//I7VfZRcB0zAOz/YDym+nt/Wka4z5R7cjNReQbNWqUrlq1qtc+76VVRUw6Iof+o7exaUV/q7s2O/V9YQ17DRcScmopXJtFaqIlEOM+EflYVUe1LO9Ug7mITBSRP4vITBEZHCw7RkT6d1eg0WLiyFyyRhSx9dN8Pt5Q4XY4JgLc/mgJTRUpHPOrGkscJux1dCXBPiLyIfAScCkwDWgeuXQZ8OeeCS9yxXk9TDijgUBdArfPteRh2hcIKM/Mj0MSmrh+Rqrb4RizUx2987gHSANGBLfQOpg3geO7Oa6ocP1lWXhSGnjthUR8fpuuxLRtyZoSir4oIP+AQsbta72sTPjraPKYAPw/VV3PL3tgtbtkbCwbOTSD/gcXUf51Pq9+3OpCjMYAcNvMGrQxjrMnN+HxWPuYCX+dafNoaqM8D6jrhliijogw+XwFv5e7ZtspMq2raWhi2SupxGXVcN3FdtdhIkNHk8dy4GoRCR3L0XwHMhV4u1ujiiJXT84lLqeaDxanUVlv05WYX3rkzSJqN+Qy/MgSdi+w9g4TGTqaPP4NOBRn4OB/4SSOy0XkHeAI4P/1THiRb1BOCiOOKqXux1wefb3lQozGwP2zfIBwxXQbZ2siR0fX8/gSOARnlt0pgB84E6e94zBV/aanAowGMy5zLgoPzm2r5s/Eqh9Lalm7LIfkwSVMmWDrdpjI0eE2D1X9TlUvVtX+qpqgqn1V9UJV/bYnA4wGl04oIGlQKWuX57C5zNo+zA53PFZCU1kaR55cTUaSrRhoIofNqtsLMpLiOWx8FU2lady7wHpdGYeq8tQTHiTexnaYyNPerLqdagRX1eO6Hk70+v30ZN552M8Tjwt/m6E2XYlh+ddlbF9dQO5+RZxwYF+3wzGmU9q78yhpsQ0HjgZSgOrg41HAMMC+Tu/EKaPzyNqrmC0f5/PZj5Vuh2PCwG0zq9DGeM4414fXxnaYCNPeSoLnNG/Aa0ApsIeqHq6qp6nq4cCeQBnwRu+EG7nivR5OmlRPoDaR2x+y6UpiXV2jnyUvpeLNqOWGy3LcDseYTutom8efgL+o6o+hhcHnNwP/0c1xRaVrp2bgSWrklecSaLLpSmLaE0sKqf4+lz3GlDCsT5rb4RjTaR1NHn2BxDZeSwCsj2EHHLpHFn1HFlG+Np/XV5e4HY5x0X2zfKDCjKnWZ8VEpo7+z10K/F1E/mVOdxE5FPg78E43xxWVRIRzzw+gTV7unF3rdjjGJVvK6vhyaTZJA0uZ/iv73mUiU0eTxwycNo8PReQnEVktIj8BHwTLZ/RUgNHmDxfkEpdVw3uvpVLdYIMGY9Fd84vxlaRz+ElVZKbY2A4TmTo6wnyzqh4MnAo8CLwffDxVVUeq6uYejDGqDM5NYdiRpdT+kMsTb9l0JbFGVZn/mAeJ83Ptr21texO5OrVcmaouAhb1UCwxY9oUD9e/Itw/18eMk92OxvSmD9eXs/XTfLL3KWLiSKuyMpGrzTsPEUkJ/XlnW++EGx2mn1xA0oAy1izLYmt5vdvhmF50y8xKAvUJTDqnkTivNZabyNXe/94qERkd/LkaqNrJZjooMyWeUSdU4ivK4L5nrOoqVtT7/Lz5YjLetDpumJrtdjjGdEl71VZTge9Cfm65gqDpgqsvT2LFowEee1z47+luR2N6w4JlRVStL2DPE7awd/9BbodjTJe0lzyGsGNsx9vAVlW11Yy6yWmH5ZExrJjNK/P4cnMl+w3McDsk08PumdkI6mHqZVZdZSJfe/+Lb2LH2uQbgJE9H07sSIzzcsJpdfhrkrh9Xrnb4Zgetq2ins+WZJHYv4xfn5bvdjjGdFl7yaMI2Cf4sxBG1VYicrSIPCAis0XkPbfj2VXXT0/Hk+jjpWfj8QfC5vSaHnDPgiJ8RRkcOr6S7NQEt8MxpsvaSx7PAg+JSDFO4lgsIoVtbR39QBGZG3zPly3KJ4jIOhFZLyJ/bO8YqrpcVX8NvAw83NHPDjeHD8um4KBCStfk8+bnNjFxtFJVHn1UwOvn2iuS3Q7HmG7RXpvHVThtHXsDf8VJJt0xGHAecC/wSHOBiHiB+4Dxwc9YKSILAS/wtxbvn6qqzcnqAmBaN8TkChHh7HMD3PthHHfNqeWke9yOyPSETzZUsGVVAVkjijhltI3tMNGhzeShqoqTMBCR44HbVfXrrn6gqi4TkaEtikcD61X1++DnPQlMUtW/Aae0dhwRGQxUqGqb3YRFZAbBqVMGDx7c1dB7xDUX5XD/X2tZviiV2tubSEno1LhNEwH+PquCQF0Wp5zdSLyN7TBRoqPTk4zrjsTRjgHAppDnm9nRWN+WacBD7e2gqjNVdZSqjsrPD89Gyj0KUtnjiBJqNuTy5FIb8xFtGpr8LH4hGW9qPTdMzXI7HGO6TYe/5opIf5y7gIFAUsvXVfXGboxrp1T1pt78vJ40bYrwb68J/5zjY+qJbkdjutPz7xZT+U0+ux37EwcMHuh2OMZ0mw4lDxE5A5iP0wZRCDS22EWBriSPLUDoqKmBwbKYcPkpBfylXzlfLM2ksKqegvRf5GYToe6a2QABD1OmuB2JMd2roxWw/wu8DvRR1QGquluLbfcuxrESGCYiu4lIAjAZWNjFY0aM7NQERh5XSWNhJv981qquokVRVQOfvpVJYt8KfntmeFabGrOrOpo8BgH/UNXSrn6giMzHmdJ9LxHZLCLTVLUJp3fXYmAtsEBV13T1syLJ1ZcngifAo4+5HYnpLvc9W0Tj9kxGHl9BXlpbC3EaE5k62ubxHrAX8GZXP1BVz2+jPKanez9jTB7pexaz8cM81v5Uxd79090OyXTRvIcBT4BrrrBBgSb6dPTO41pghohcKiL9bUr27pcU7+W4U+rwVydz+8Nlbodjumj1xgo2rcwjc68izjjCxnaY6NPR5PE5sD9O19hN2JTsPeIP09KQBB8Ln4kjYNOVRLRb5pQTqEliwhn1JMTZ2A4TfTpabWVTsveCo0fkkH/AVoo+L2DpVyUct1+e2yGZXdDYFODV55LwpDRw4/Qst8Mxpkd0KHmo6rwejsMAHo9w5jlNPLAqjjtm1XDc3ZY8ItHCj4oo/zqPwUduZeTQnY11NSYy2f10mLn20hy8GbW8syiFuka/2+GYXXDXg/Xg93LJJYqIuB2OMT2izTsPEfkImKKqX4nISnZSbaWqo9t73XTMsD5p7H7YJr59cwBPryjkkuP6uh2S6YTSmkZWvZFJQkEFvzvHxnaY6NXenccaoC7k551tpptccomAerhvdsuB/Cbc3f98IQ1bszjwuHIKMmymABO92ptV97KQn6f0SjQGgF+fns9/9a1g9ZJMiqsbbIBZBJk7T8ET4OoZ9m9mopu1eYShvLREDhpXQeO2TB580aYriRRrtlSy8cN80vcs5uyjrMrKRDdLHmHqyukJIAHmPWw9pCPFLXPK8FcnMf70OpLivW6HY0yPsuQRps49Op+03UvY8GEu32yrdjscsxNN/gAvPZOIJ6mRG6dnuB2OMT3OkkeYSor3cuyvavFXpnDno12ej9L0sEUfF1P2VT79Dilk9J5ZbodjTI+z5BHGrr08FUlo4rkFNl1JuLvjwTrwe7nwYhvbYWKDJY8wNnafXPL2K6To83yWf213H+GqvLaRDxdnEJ9Xxe/Ps1kBTGyw5BHGPB5h0tlNaGM8d86xdo9wNfPlIuq3ZLPf2DL6ZyW7HY4xvcKSR5i79tJsvGn1vP1yMvU+m64kHM2eEwBRfjcj3u1QjOk1ljzC3N790xkyuoiq9Xk8/16x2+GYFtZtreL79/NI26OYycfauh0mdljyiAAXXwwEPNwzu97tUEwLtz5Uhr8qmeNOrSU5wcZ2mNhhySMC/PasfBIKKvjkrUxKa2y+q3DhDygvPB2PJPq48QpbNtjEFkseEaAgPYn9jq2gYWsWs18udDscE7R4dQmlXxbQ96BCxgzPdjscY3qVJY8IceX0eBBlrk1XEjbueLAGbfIy+cKAje0wMceSR4SYPDaf1N1K+O69XL4vqnE7nJhXWe/jvVfTicup5g8X5rodjjG9zpJHhEhJiOOoiTU0VaRwx2MlbocT8+YuKqJuUw77HFvKoJwUt8MxptdZ8ogg105PReKbePZJL6pWfeWmB+f4AXWqE42JQZY8IsjxB+SSs08RhZ8V8MG3ZW6HE7O+K6zh23dzSd29hItOsHU7TGyy5BFBvB7h1LN8BBriuW1OldvhxKzb5pXgr0jh2F/VkJLQ5mKcxkQ1Sx4R5topWXhT63lzYRINTTZdSW8LBJRnn4pDEnzccEWa2+EY4xpLHhFm/0EZDDy0mMpv8nnxA5uupLe99UUJxV8UkH9AEcfsneN2OMa4xpJHBLrwQoWAh3vn2HQlve32mTWoL45zL2jC47GxHSZ2WfKIQFednU98XhUr38ygotbndjgxo7qhieWvpBGXVcP1F9u6HSa2WfKIQP2yktj3mDLqN2cz51WbrqS3PPx6IbUbc9nr6BKG5NnYDhPbLHlEqN9MiweU2Q9Zo3lvuX9WEwC/nmY9rIyx5BGhLjw+n5ShpXz7bi4bi2vdDifq/VhSy7oVuaQMKWHKBFu3w5iITB4iso+ILBCR+0XkbLfjcUNqYhxjTqqmqTyVu+fbdCU97bZHi2kqS+Wok6tJS7Q7D2N6PXmIyFwRKRSRL1uUTxCRdSKyXkT+uJPDTATuUdXfAJf0WLBh7prLU5A4P0/N99h0JT0oEFAWPBGHxDdx3YxUt8MxJiy4cecxD5gQWiAiXuA+nKSwD3B+8O5ifxF5ucVWADwKTBaRW4GYndL0pINyydq7iG2f5rPyu3K3w4lay9eWUvhZPnn7F3L8ATH7382Yf9HryUNVlwGlLYpHA+tV9XtVbQSeBCap6heqekqLrTC4XQn8EWhzpJyIzBCRVSKyqqioqMd+J7fEeT386sxGAvUJ3D630u1wotatM6vRxnjOPK8Jr43tMAYInzaPAcCmkOebg2WtEpGhIjITeAS4ta39VHWmqo5S1VH5+dE5gd21UzLxpDSw+MVEfP6A2+FEndrGJpa+nIo3o5brL7UR5cY0C5fk0Smq+oOqzlDVC1V1hdvxuOmgIRkMOKSIinX5vLLSpivpbo+/VUTN97kMO7KYPfvYXFbGNAuX5LEFGBTyfGCwzOyEiHD+hQp+L3fPrnM7nKhz32wfIMywsR3G/ItwSR4rgWEispuIJACTgYUuxxQxrj43j/jcaj54PZ2KOpuupLtsLqtjzdIckgaVMvXk6Kz2NGZXudFVdz7wPrCXiGwWkWmq2gRcBSwG1gILVHVNb8cWqQZkJzPiqFLqN+XwyOs2XUl3ueuJIppK0xgzoYrMZFsx0JhQbvS2Ol9V+6lqvKoOVNU5wfJFqjpcVfdQ1f/p7bgi3YzLnGqVmXNtupLuoKrMf8yLxPm59gqbx8qYlsKl2sp00aUn5ZM8uJSvV+SwqdSmK+mq978pY9unBeTsW8iEkTaDrjEtWfKIEulJ8Rx+YhVNpWnc85RNV9JVt8ysItAQz6RzG21shzGtsOQRRa65PBm8fuY/ITZdSRfU+/y8tTAFb3odN0yxsR3GtMaSRxSZeEgeWSOK+emTfD79wUac76r5Swup/i6P3Y8oZkT/dLfDMSYsWfKIIvFeDxMmNRCoTeT2h8rdDidi3TvLBypMv8yqq4xpiyWPKHPd1Aw8yY28+nwiTTZdSadtLa/niyXZJA0oY8Zptm6HMW2x5BFlDtk9k34jiyj/Op9XP7GG8876x1NF+IrTGX1iJVkpCW6HY0zYsuQRZUSEyRf40SYvd8+2Lrudoao89qiA18+1VyS7HY4xYc2SRxS6enIecdnVvL84jeqGJrfDiRgrv6vgp48LyN67iJNH2dgOY9pjySMKDc5NYa8jy6jdmMOjb9h0JR11y6wKAvUJnHJWI/Fe+9Mwpj32FxKlpk/1AMKDc+3OoyPqfX7eeDEZb1o9N0zLcjscY8KeJY8oddmEApIGlvLV8mx+Krep2nfm2RXFVH6bx5DRRew/KMPtcIwJe5Y8olRmcjyHnlCFrzid+56xRaJ25h+zGiDg4bIpNrbDmI6w5BHFrr48Cbx+Hn/MpitpT2FVPavfyiKxbzm/OcPW7TCmIyx5RLHTRueTObyYTavy+HyTTVfSlnsXFNNYmMEh4yvITUt0OxxjIoIljyiWEOdh/Gn1BGqSuMOmK2mVqvLwI4AnwDVXJLkdjjERw5JHlLtuWgaepEZefi4Bf8CqrlpavbGSzavyyRpRxKTDrMrKmI6y5BHlDtsziz4HFlH2VQFvrLbpSlq6ZVYFgdpEJpxRT0Kc/TkY01H21xLlRISzJzvTldw5p8btcMJKY1OAV59PxJPSwI3Ts9wOx5iIYskjBlx7US5xWTWseDWVGpuu5GcvfFBExbp8Bh1axEFDbGyHMZ1hySMGDM1LZc8xpdRuyOWJJTZdSbO7HqyHgIcpU5w7NGNMx8W5HYDpHVOneLhxkfDAXB+XT3A7mo5r8geo8/mpa/RT2+inzuc8Vtf5KasMUFbhp7JaqagMUFEFNTVKVbVSXQO1NVBbC3V1Ql0tNNQJDfUeGusFX4OH6h/7kNCngivPsoZyYzrLkkeMmH5yPn/uX8aXS7MprKynIKN7uqX6/AFqG/3UBy/qdY1+6nxNVNb4Ka8MUF4ZoKzSuaBXVgUv7NVKdW3Ixb1WqK+D+jrnwt4QvLg3NXjwN3rRJi8Bnxdt9KJN8QR8XvB7OxmpIvH+nzdPvB9vRh2nXFJBfvrQbjkXxsQSSx4xIjs1gUOO38p7jw7hrw9/xxnjMiirDFBRFaC8UqmsVqqqAlRWORf2mmrnwr7j4u5szoVdaKz30NToIdAYvLD7vKgvDvUlOBf3QGdrRBVJ8CPxTXji/Ujcjot8fEoT8YkB4pMCJCYpiclKUpKSnKIkp0BqCqSmOltGmpCWJmRkCFnpHjLThexMD5lpXlITvaQkeEmK95KSkEhyfArJCdk9cr6NiXaWPGLI1Zcn8t7jAe67enfuo5N1/BL4+WLe8ht8UqpzcU9ICpDQ4uKekgopyZCWBqmpQnqac4FPTxOyMjxkZTiPmWlxpCR4SU5wLvDJ8fEkJ3hJjPNYe4QxYUhiZc4jESkCNrocRh5gsxQ67FzsYOdiBzsXO4TLuRiiqr9oGIyZ5BEORGSVqo5yO45wYOdiBzsXO9i52CHcz4V11TXGGNNpljyMMcZ0miWP3jXT7QDCiJ2LHexc7GDnYoewPhfW5mGMMabT7M7DGGNMp1nyMMYY02mWPHqBiJwjImtEJCAio1q89u8isl5E1onISW7F2JtEZELw910vIn90O57eJCJzRaRQRL4MKcsRkTdE5NvgY0wMexeRQSKyRES+Cv59/D5YHnPnQ0SSROQjEfkseC7+M1i+m4h8GPxbeUpEEtyOtZklj97xJXAmsCy0UET2ASYD+wITgH+KSGcnbYoowd/vPmAisA9wfvA8xIp5OP/Wof4IvKWqw4C3gs9jQRNwnaruAxwOXBn8vxCL56MBOE5VDwQOAiaIyOHA34E7VXVPoAyY5l6I/8qSRy9Q1bWquq6VlyYBT6pqg6puANYDo3s3ul43Glivqt+raiPwJM55iAmqugwobVE8CXg4+PPDwOm9GZNbVHWrqn4S/LkKWAsMIAbPhzqqg0/jg5sCxwHPBMvD6lxY8nDXAGBTyPPNwbJoFou/8870UdWtwZ+3AX3cDMYNIjIUGAl8SIyeDxHxishqoBB4A/gOKFfV5hXcwupvxSZG7CYi8ibQt5WX/qSqL/Z2PCYyqaqKSEz1nxeRNOBZ4BpVrQydCDOWzoeq+oGDRCQLeB4Y4W5E7bPk0U1U9YRdeNsWYFDI84HBsmgWi7/zzmwXkX6qulVE+uF884wJIhKPkzgeV9XngsUxez4AVLVcRJYARwBZIhIXvPsIq78Vq7Zy10JgsogkishuwDDgI5dj6mkrgWHBXiQJOB0GFrock9sWApcGf74UiIk7VXFuMeYAa1X1jpCXYu58iEh+8I4DEUkGxuO0AS0Bzg7uFlbnwkaY9wIROQO4B8gHyoHVqnpS8LU/AVNxep5co6qvuhVnbxGRk4G7AC8wV1X/x92Ieo+IzAfG4ky3vR24CXgBWAAMxlk24FxVbdmoHnVE5ChgOfAFEAgW/wdOu0dMnQ8ROQCnQdyL86V+gar+VUR2x+lUkgN8Clykqg3uRbqDJQ9jjDGdZtVWxhhjOs2ShzHGmE6z5GGMMabTLHkYY4zpNEsexhhjOs2Sh4loIjJFRDQ4Srkz77tZRIp7MK4fROS2Hjz+KhGZ11PHN2ZnbIS5MT3jDKDE7SCM6SmWPIzpAar6qdsxdJfgFCKB4NxLxgBWbWXCmIgcISILRWSriNSIyGoRuXAn7xkarMa6QEQeFZGq4OJLN7Wx/0gR+UBEakXkUxE5usXrl4jIChEpFZGy4OJFo1o7Vov3/Uu1lYjMC1Y1jReRz4O/zwoR2bcDx9pPRN4VkXoRWSsip7Wx39Ei8k7wdykRkVkikt5in7HBz68XkZUiMlpEikXk5pB9lorIMyIyQ0S+A+qB/sHXpgcXK2oQkY0icuOuxGEin915mHA2BHgXeADnAnYk8JCIBFR1/k7eeyvwMs68QMcAN4lIsareF7JPCs6UEHfiTP19E/CciAxR1drgPkOBR3Cmx04AzgeWi8i+qvp9J3+fwcG4/geoA24DnhKR/bWNqR6C8xwtBoqBC4BknKld0nAWGWve70jgTZypTs4GcoH/A7KDzxGRAcAi4D2caUD6Ao8Hj9nSkcAewL8BtUCFiNwA/C9wC7AUOAT4LxGpVdV7OxqHiRKqapttYb8BgvNl50Hg7ZDyKTiL5qQFnw8NPn+9xftn4cxI6gk+vzm433Eh+xwULJvQRgyeYAxfA3/ZSbw/ALeFPJ+HM3/ZsJCy04OfN6Kd4/wW8AEDQ8qODL5vXkjZcmBJi/ceF9xvv+DzW3GSUHLIPucG97k5pGwpTnLrE1KWAVQDN7X4jL/iJF5vR+OwLTo2q7YyYUtEskXkHyKyEecC6gNmAMM78PbnWzx/DqfqZWBIWSPOhbLZV8HHn/cRkb1F5HkR2Q74gzHs1cEYWvpBVb9t7/NaMRr4WFU3Nxeo6ruETFMuIik403cvEJG45g1YEYz3kOCuhwJvqGpdyPHbmtH4Y1XdHvL8CCAVeLrFZ7yNs1jTwE7EYaKAJQ8TzuYB5+F8Yz4R5+I3F0jqwHtbrgHR/LxfSFmVqjbP5oo6y+LSfPxgPf3rOOuPXAscHYzhsw7G0FJ5i+f/8nlt6Evr61mElmXjzMb6T3YkWR/Outjx7Fg/pS9QFHoQVa3HuaNoaXuL53nBxzUtPmNJsHxQJ+IwUcDaPExYEpEk4BTgSlV9IKS8o194Ctp4vrXlju04AueuYLyqfh0SQ2YnjtFV22h9RbnQ36+cYNUTTptGSz+FHCs/9IXgeW5tjEzLNpjmKdFP4ZeJBWAdzrTqHYnDRAFLHiZcJeLcGf+8dkHwTuA0fnlha80ZwP0hz8/ESRybW9+9Vc0NyaExjMFpV/m4E8fpipXAhSIysLnqKtgo/XPyUNUaEfkA2EtV/7qTY10mIskhVVet9txqxfs47SD9VfWVtnbqYBwmCljyMGFJVStEZCXwFxGpxPlW+0egAqfxdmf2FZEHcZY4PQaYBvw+tJqqAz7AqdKZJSK34NyF3EzvLgX6EPD/gFeC3WmTgf/CafgOdSPwlogEgGeAKpzeXb8C/qSq3+D00roSeElE7sSpxvojTm+qds+LOkuj3gzcLSJDgGU4yX04ME5Vz+hEHCYKWJuHCWcXAN/jdJW9GycRPNLB996Ik2SeBa7AueDe25kPDzYYn4NzkX0RuAb4NbC+M8fpCnW6DJ8E1OCsKHcTcB3OCnuh+63ASZL5wKPASzjnYBPBaiZV3YJzES/A6UDwO5xVLL1AZQdiuQWnw8JEnPMxH7gQp4dVh+Mw0cFWEjRRRUSGAhuAU1X1ZZfDCXuyYynY41R1yc72N6aZVVsZE0NE5O84a2Fvw+ly/Gfgc+AdN+MykceShzGxJRGn63MfnPaI14FrO9kWZIxVWxljjOk8azA3xhjTaZY8jDHGdJolD2OMMZ1mycMYY0ynWfIwxhjTaf8f2BiMvazq4XEAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"#data\n",
"x = np.linspace(-15, 35, len(Laverage))\n",
"\n",
"fig, ax = plt.subplots()\n",
"\n",
"ax.fill_between(x, Lmin, Lmax,alpha=.5,linewidth=0)\n",
"ax.plot(x,Laverage,linewidth=2)\n",
"\n",
"ax.set(xlim=(-15, 35),ylim=(min(Lmin), max(Lmax)))\n",
"ax.set_yscale('log')\n",
"\n",
"\n",
"plt.plot(x,Laverage,'b')\n",
"plt.title('BS-based', fontsize=20)\n",
"plt.xlabel('alpha in degree', fontsize=15)\n",
"plt.ylabel('infidelity 1-F', fontsize=15)"
]
},
{
"cell_type": "markdown",
"id": "e35027ee",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"## 2. Perceval implementation of the MZI-based interferometer\n"
]
},
{
"cell_type": "markdown",
"id": "d3586f39",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"Now we do the same construction for the MZI-based circuit. The main difference between both circuits is the structure of the building blocks. Note that the MZI-based circuit has $m$ building blocks with four components whereas the BS-based circuit has $2*m$ building blocks with two components."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1d6abe25",
"metadata": {
"pycharm": {
"is_executing": true,
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"def mzi(P1,P2):\n",
" mz= (pcvl.Circuit(2, name=\"mzi\")\n",
" .add((0, 1), phys.BS(theta=(45)/180*np.pi))\n",
" .add(0, phys.PS(P1))\n",
" .add((0, 1), phys.BS(theta=(45)/180*np.pi))\n",
" .add(0, phys.PS(P2)))\n",
" return mz"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e7fd4aa8",
"metadata": {
"pycharm": {
"is_executing": true,
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"c = pcvl.Circuit.generic_interferometer(N,\n",
" fun_gen=lambda idx: mzi(pcvl.P(\"phi_m%d\" % idx),pcvl.P(\"phi_n%d\" % idx)),\n",
" shape=\"rectangle\",\n",
" depth=N,\n",
" phase_shifter_fun_gen=lambda idx: phys.PS(phi=pcvl.P(\"phi_r%d\" % idx)))\n",
"pcvl.pdisplay(c)"
]
},
{
"cell_type": "markdown",
"id": "f08b44d5",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"### Basinhopping algorithm"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c0afc0a7",
"metadata": {
"pycharm": {
"is_executing": true,
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"def mzi(P1,P2,alpha):\n",
" mz= (pcvl.Circuit(2, name=\"mzi\")\n",
" .add((0, 1), phys.BS(theta=(45+alpha)/180*np.pi))\n",
" .add(0, phys.PS(P1))\n",
" .add((0, 1), phys.BS(theta=(45+alpha)/180*np.pi))\n",
" .add(0, phys.PS(P2)))\n",
" return mz"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ff2dc28f",
"metadata": {
"pycharm": {
"is_executing": true,
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"def calculate_angle(index_unitary, U_target, alpha):\n",
" start = time.time()\n",
" c = pcvl.Circuit.generic_interferometer(N,\n",
" fun_gen=lambda idx: mzi(pcvl.P(\"phi_m%d\" % idx),pcvl.P(\"phi_n%d\" % idx),alpha),\n",
" shape=\"rectangle\",\n",
" depth=N,\n",
" phase_shifter_fun_gen=lambda idx: phys.PS(phi=pcvl.P(\"phi_r%d\" % idx)))\n",
" params = c.get_parameters()\n",
"\n",
" # goal is to find phis matching U0 - by minimizing fidelity\n",
" infidelities = []\n",
"\n",
" for _ in range(n_try):\n",
" init_params = np.random.randn(len(params))\n",
" res = basinhopping(lambda x: infidelity(c, U_target, params, x), init_params, stepsize=0.1, niter=n_iter)\n",
" infidelities.append(res.fun)\n",
" \n",
" return min(infidelities)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0f75ba07",
"metadata": {
"pycharm": {
"is_executing": true,
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"n_angles = (abs(angle_min)+abs(angle_max))/angle_step\n",
"\n",
"def discovery_unitary(index_unitary):\n",
" np.random.seed(int(10000*time.time()+os.getpid()) & 0xffffffff)\n",
" U_target = pcvl.Matrix.random_unitary(N)\n",
" l_infidelities = []\n",
" M=angle_min + math.ceil((angle_max-angle_min)/angle_step)*angle_step + 1\n",
" for alpha in tqdm_notebook(range(angle_min,M,angle_step), desc = \"Angle progress bar unitary {}\".format(index_unitary+1),leave=False):\n",
" l_infidelities.append(calculate_angle(index_unitary, U_target, alpha))\n",
" return l_infidelities\n",
"\n",
"try:\n",
" os.remove(\"%s-%d.log\" % (logfilemzi,N))\n",
"except OSError:\n",
" pass\n",
"\n",
"for index_unitary in tqdm_notebook(range(n_unitary), desc = \"Unitary progress bar\"):\n",
" result = discovery_unitary(index_unitary)\n",
" e = datetime.datetime.now()\n",
" with open(\"%s-%d.log\" % (logfilemzi,N), \"a\") as f:\n",
" f.write(e.strftime(\"%Y-%m-%d %H:%M:%S\")+\"\\t\") \n",
" f.write(\"\\t\".join([\"%g\"%inf for inf in result]))\n",
" f.write(\"\\n\")"
]
},
{
"cell_type": "markdown",
"id": "0a40d8ba",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"## Result of the MZI-based scheme"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "38fa883b",
"metadata": {
"pycharm": {
"is_executing": true,
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"with open(\"mzibasednotebook-opt-5.log\", \"r\") as filin:\n",
"... for ligne in filin:\n",
"... print(ligne)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d5926a3c",
"metadata": {
"pycharm": {
"is_executing": true,
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"L=[]\n",
"with open(\"mzibasednotebook-opt-5.log\", \"r\") as flog:\n",
" for l in flog:\n",
" L.append([float (f) for f in l.strip().split(\"\\t\")[1:]])\n",
"\n",
"A=np.asarray(L)\n",
"print(A)\n",
"print(A.shape)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fe4bdf97",
"metadata": {
"pycharm": {
"is_executing": true,
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"Lmin=[]\n",
"Lmax=[]\n",
"Laverage=[]\n",
"\n",
"for i in range(A.shape[1]):\n",
" Langle = A[:,i]\n",
" Laverage.append(np.average(Langle))\n",
" Langle.sort()\n",
" Lmin.append(np.average(Langle[0:10]))\n",
" Lmax.append(np.average(Langle[-10:]))\n",
"\n",
"print(Laverage)\n",
"print(Lmin)\n",
"print(Lmax)"
]
},
{
"cell_type": "markdown",
"id": "c8c6dd9c",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"### Graph of the infidelity as a function $\\alpha$ in degree : MZI-based"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d4502f9e",
"metadata": {
"pycharm": {
"is_executing": true,
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"#data\n",
"x = np.linspace(-15, 35, len(Laverage))\n",
"\n",
"fig, ax = plt.subplots()\n",
"\n",
"ax.fill_between(x, Lmin, Lmax,alpha=.5,linewidth=0)\n",
"ax.plot(x,Laverage,linewidth=2)\n",
"\n",
"ax.set(xlim=(-15, 35),ylim=(min(Lmin), max(Lmax)))\n",
"ax.set_yscale('log')\n",
"\n",
"\n",
"plt.plot(x,Laverage,'r')\n",
"plt.title('MZI-based', fontsize=20)\n",
"plt.xlabel('alpha in degree', fontsize=15)\n",
"plt.ylabel('infidelity 1-F',fontsize=15)\n"
]
},
{
"cell_type": "markdown",
"id": "90dd9c90",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"## Our results"
]
},
{
"cell_type": "markdown",
"id": "e6c23a2d",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"For matrices of size 5, we obtained this graph below."
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "b9f496e7",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
""
]
},
{
"cell_type": "markdown",
"id": "2457fba0",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"The MZI-based interferometers are equally sensitive to both positive and negative values of α. The range of tolerance for a correct implementation is of several degrees. For the BS-based interferometers, infidelity behaves rather differently: while for α < 0 the performance of the two are comparable, for α > 0 the BS-based interferometers provide a correct implementation for a much larger range of α (as\n",
"large as ∼ 20 degrees). Therefore, the BS-based circuit is more stable than the MZI-based circuit although the BS-based implementation fail at implementing some unitary matrices (like permutation) in the ideal case."
]
},
{
"cell_type": "markdown",
"id": "3130c4e7",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"In the article [1], they also performed optimisation on matrices of size 10. We pushed the calculations a bit further on matrices of size 12."
]
},
{
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
""
]
},
{
"cell_type": "markdown",
"id": "1c464c9d",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"We only made our computations for the BS-based scheme because we wanted to analyze the stability of this new circuit. We do not have done both yet because it is very time consuming. Indeed, for only one matrix of size 12, the computation takes about 37 hours. Since in our code, we loop over 300 matrices, we get a total time it would take to complete this calculation of roughly 11,100 hours (~460 days) of calculation in total."
]
},
{
"cell_type": "markdown",
"id": "a516779b",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"To conclude, the BS-based circuit is a very stable scheme but not universal. MZI-based circuit is universal when BS are perfectly manufactured, although it is not possible in practice. However, as one can see on the graphs, the MZI-based circuit is much less robust to errors. In practice, it is common to have $\\alpha$ in the range of 5 degress. Thus, the MZI-based circuit is not really universal in practice.\n",
"\n",
"It is therefore necessary to continue to seek relevant criteria to properly evaluate the interferometer schemes."
]
},
{
"cell_type": "markdown",
"id": "c00d6e3b",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"## References"
]
},
{
"cell_type": "markdown",
"source": [
"[1] Suren A Fldzhyan, M Yu Saygin, and Sergei P Kulik. Optimal design of error-tolerant reprogrammable multiport interferometers. [Optics Letters](https://doi.org/10.1364/OL.385433), 45(9):2632–2635, 2020.\n",
"\n",
"[2]: Michael Reck, Anton Zeilinger, Herbert J Bernstein, and Philip Bertani. Experimental realization of any discrete unitary operator. [Physical review letters](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.73.58), 73(1):58, 1994."
],
"metadata": {
"collapsed": false,
"pycharm": {
"name": "#%% md\n"
}
}
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.4"
}
},
"nbformat": 4,
"nbformat_minor": 5
}